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Abstract 

This paper presents generic models for the effect of a chemical 
toxin on cells forming the tissue of an organ. The models are 
illustrative, not specific to organ or toxin. Interactive response of 
within-tissue toxin and cells is modeled: cell capability to modify 
(metabolize or bind) toxin is represented, as is the alteration of that 
capability by toxin presence. Both processes are represented in the 
context of simple versions of the cell cycle. 

Toxin-cell interaction is explicitly represented in terms of inter- 
cell chemical signaling that encourages replacement or repair of 
cells. The impact of cell death by apoptosis and necrosis is 
represented in terms of within-tissue toxin concentration. The 
model is explicitly stochastic, representing inter-cell and toxin input 
and within-organ concentration in terms of diffusion 
approximations. Explicit mathematical discussion is given of dose- 
response function behavior at low doses. 



1. Introduction 

Application of laboratory toxicology data to environmental and human 
problems of risk assessment almost always requires extrapolation of the data 
from the experimentally-used dose regimen to the exposure conditions of 
practical concern, and from the animal species tested to the species of concern 
(usually man). This extrapolation and estimation process is known as chemical 
risk assessment. The risk assessment process has undergone considerable 
evolution, moving from a qualitative basis for decision making to an increasingly 
quantitative basis, and from the use of default assumptions to the application of 
mechanistic mathematical models as tools upon which to base decisions. In the 
context of determining safe human exposure limits to potentially toxic chemicals, 
there are two sub-tasks to be accomplished: estimation of low-dose risk in 
animals, and the subsequent conversion of animal risk estimates to human risk 
estimates. An authoritative survey of many of the statistical issues and 
opportunities is given by Krewski and Franklin (1991). 

Extrapolation to low-dose effect, the first sub-task, cam be accomplished either 
by assuming that biological response varies in a specified mathematical manner 
(e.g. probit, logit) with organ host exposure, or by using physiologically-based 
pharmacokinetic (PBPK) models to relate organ dose to host exposure; cf. 
Krewski and Franklin (1991), chapters 8 and 9; M. Andersen, H. Clewell Jr., and 
C. Travis have been prominent in this research area. In the latter approach, multi- 
compartment physiological models are formulated using actual tissue volumes 
from the experimental and target species and actual perfusion rates to provide 
for chemical transport between the compartments. Thus the pharmacokinetics of 
high to low-dose extrapolation become amenable to calculation, and external 
measures of external dose or exposure Cem be translated to concentrations in the 
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target organ (internal dose). Target-organ chemical concentration may be 
translated into estimates of risk if a suitable biologically acceptable mathematical 
model can be used to relate chemical presence in the organ to harmful outcome. 
Such models, called pharmacodynamic, exist for carcinogens, in the form of the 
widely-used linearized multistage approximation for cancer dose/resp>onse, and 
multistage clonal expansion models; cf. Moolgavkar (1988). These models are 
based on the general concept that chemical alteration of the cellular genes may 
give rise to permanent, heritable changes in the genetic information stored in the 
cell nucleus, and lead to phenotypic changes in the altered cells that ultimately 
cause the formation of malignant tumors. Simulation is sometimes used to study 
stochastic models of carcinogenesis for large numbers of cells; (cf. Bois (1992)). 

Unfortunately, analogous dose/response models are not widely available for 
toxic responses other than carcinogenesis. The mathematical expression of even 
the relatively simple concept of genetic alteration leading to cancer involves 
significant simplification of biological reality, and significant mathematical 
complexity. This state of affairs is exacerbated when one attempts to describe the 
interactions leading to loss of cell function and cell death in tissues of a whole 
animal. The multi-layered control and response systems present in an intact 
living animal are poorly understood and thus have not yet been adequately 
modeled. Nonetheless, these control mechanisms defend against the majority of 
toxic effects observed as a result of chemical exposure, and their failure 
represents much of what is observed as the expression of toxicity. 

This paper describes initial formulations of models for cellular response to 
toxic insult. Beginning with a simple Markov model, successive models increase 
in biological resolution or refinement as well as mathematical complexity. Our 
model formulation takes into account broad aspects of the current state of 
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knowledge concerning the control eind regulation of cellular properties in tissue. 
The design of these models is not to provide a detailed description of the 
response of specific organs to toxic insult, but rather to begin the process of 
mathematically describing recent significant advances in knowledge that have 
been made in cell biology, hormonal regulation, and hormetic control mech- 
anisms that regulate cell response to toxicants. Some modifications of the model 
formulations may be needed to apply them to specific organs; for one thing the 
complex geometry of organ architecture is not represented. To have included 
complex, three-dimensional relationships between different cell types would 
have increased the mathematical complexity considerably. But even in its present 
rudimentary state the model is useful for stimulating discussion and suggesting 
experimentation, with the goal of furthering the understanding of cellular 
response to toxic substances including the effects of tissue-mediated responses 
(e.g. mitogenic stimulation). The current formulation of the model, or extensions 
thereof, is suitable for experimental evaluation in cell culture if the experimental 
conditions are properly chosen and chemical-spedfic details are added. 

We provide the following abbreviated biological background. The tissue 
making up organs almost always involves a complex geometrical juxtaposition of 
several cell types having specific and often overlapping functions. This tissue 
architecture may be maintained by the presence of a nonliving matrix of proteins 
collectively known as a basement membrane. Interactions between this membrane 
and the cells are known to be critical to its stability as a mature, functioning 
entity. The whole of the tissue is permeated by branching blood vessels, each 
generation of which is successively smaller; these serve to provide a constant, 
stable milieu in which the cells exist. Nutrients, control signals in the form of 
specific biomolecules, and xenobiotic chemicals are brought to the immediate 
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surroundings of the cell by the vasculature in the tissue and cellular metabolic 
products; the products of energy metabolism and cellularly derived control 
biomolecules are removed from the cellular microenvironment by the same 
means. 

The state of the cell at any time is a reflection of its age, the summation of the 
control chemicals reaching and leaving the cell, the effects of xenobiotic 
chemicals present, if any, and the state of the cells surrounding it. Cellular 
contact with the surrounding cells and basement membrane act to supply 
chemical signals that modulate its activity. A given cell may be (i) nearly 
quiescent, (ii) active biochemically, i.e., producing metabolites of absorbed 
materials for its own use or for export, (iii) in a state of stress due to shortage or 
excess of biochemical molecules, (iv) in a process of programmed cell death 
(apoptosis), (v) in the process of dying from chemical insult (necrosis), or (vi) 
dividing to form progenitor cells in response to a need to replace cells already 
lost; these conditions need not be mutually exclusive at any point in time. Some 
specific cells may alter or completely change their observable characteristics in 
response to chemical signals. The most well-known of these cells are the 
pluripotent stem cells of the hematopoietic system. 

The models to be presented here represent some, but by no means all, of the 
features mentioned. They provide a basis for proceeding further. 

2. Model Structure 

In the following sections of this paper we present a sequence of mathematical 
models that represent interactions between the cells in an organ's tissue and a 
toxic chemical or agent (hereafter called toxin) originating externally and entering 
the organ and coming into contact with the cells, possibly in modified form. The 
initial models postulate a mature organ with capacity to contain, at most, a fixed 
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number of cells. Broadly speaking, the organ's cells attempt to remove (e.g. 
metabolize or bind) the toxin, but are, to a dose-dependent degree, also affected 
by the toxin concentration. We distinguish between cells that are in one of two 
stages of the cell cycle: these we call, oversimplifying the true cell cycle, 
functional, and dividing (S-phase and mitotic). Functional cells are assumed to be 
capable of removing toxin, but are also susceptible to premature death because of 
the toxin's action. Dividing cells are in the process of DNA replication and actual 
mitosis, here splitting into two daughter cells; while in this stage the cells are 
assumed incapable of removing toxin. However, toxin presence is assumed to 
affect the cell cycle, e.g. by shortening the functional, and lengthening the 
dividing, stages. The model proposed emphasizes inter-cellulcir signaling: when 
a cell dies, i.e. ceases to function, it effectively requests one other fimctioning cell 
to divide to produce its replacement. It is biologically plausible that signaling 
occurs only to neighboring cells, but the present model, being without spatial 
characteristics, does not literally recognize this restriction. Also, signals are 
assumed to be specific: one signal is emitted per cell death, so eventually one 
replacement occurs. Alternative formulations are sometimes plausible and are 
easily studied. 

3. Single-Stage Markov Model 

We now present several alternative simple Markov models for cell-toxin 
interaction, emphasizing the representation of signaling. Although cells can age 
and exist in different stages of functionality we omit an account of this at first, 
but return to it later. 

Let D(f) represent the number of functional cells in the organ at time t, and 
M(f) the number of dividing (mitotic) cells. If a functioning cell that dies 
immediately induces another to begin division (instant signaling) then 
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Co = D(f) + 2M(0 (3.1) 

where Co is the maximum organ size, here assumed constant; realistic 
elaborations are possible. Note that organ size is, strictly sp>eaking, D(f) + M(f), 
which is the number of cells present at time f; Co acknowledges that space for all 
cells must leave room for each cell currently in naitosis to eventually become two 
differentiated cells. We describe {D(0, M(f)} as a birth and death process in 
continuous time, with D{t)e (0, 1, 2, ..., Co). Typically Co is large, i.e. of order 
10 ^ 0 - 10^2 for mature humans, and this will be exploited to carry out an 
asymptotic analysis. 

Let T(t) denote the concentration of toxin in the organ at time t. We represent 
{T(0} as a diffusion process whose drift and diffusion coefficients are both 
influenced by D(0 so as to represent metabolic or binding action, hence toxin 
removal, and by the toxin input to the organ, zff), as well as the instantaneous 
concentration, T(t), itself. 

In the following we describe specific models, beginning with the most simple 
and (potentially) least realistic, but also the most parsimonious. We explicitly 
include the effect of inter-cellular signaling. 

3.1 Cellular Signaling 

Suppose a cell dies from some cause: apoptosis or necrosis. The effect of its 
signaling to other cells for replacement can be captured by postulating that a 
ghost (release of growth factor) of the newly-dead cell now exists, the fimction of 
which is to find and interact with a neighboring live cell, inducing that cell to 
divide. The ghost then disappears after its purpose is served. One can thus 
represent the delay in mitotic response to cell death by adjusting the ghost's 
search rate for functional cells. 
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The following describes a possible Markov generator that includes the 
signaling (ghost) influence. Let G(t) denote the number of signals present in the 
organism at time t, and introduce the constraint 



Co = D(0 + G(f) + 2M(0. (3.2) 

The following then describes the possible transitions allowed to occur in a small 
time interval, (t, t + dt). 



Signaling Model 
(Markov Generator) 



t 

Tit) 



t + dt Probability 

Dit) - 1, Git) + iMt) X*Dit)dt 

Tit) + 5* dt + CfdWit) 

(cell death, signal creation) 



Dit)-1, Git)-l,Mit) + l 
Tit) +6*dt + odWit) 

(cell signaled to divide, signal 
disappears) 



e*Git)iDit)/Co)dt 

(3.3) 



Dit) + 2, Git), Mit) - 1 /z*M(f) dt 

Tit) + S*dt + odWit) 

(two newly-divided cells emerge) 



All rates of cell transition, namely X*, 6^, and are presumed to be influenced 
by current toxin concentration, 7{t). In turn, the latter is influenced by the 
metabolic or binding capability of the cells to remove toxin; the net mean 
increase of toxin in the short run is 5*(T(f), D(f),t(f)). Our subsequent analysis 
does not, in principle, depend upon specific functional forms for any of the above 
parameters; adequate smoothness and differentiability is assumed where 
necessary. When specific functions are required we shall use these: 
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(3.4,a) 



A'(T(())=a/’'<'> 









(3.4,b) 



e\T{t)) = $oe^'^^^^) A*,0‘,/>O. (3.4,c) 

Note that it has been assumed that increased toxin level increases cell death rate, 
decreases the rate of completing mitosis, and slows the rate at which ambient 
signals reach their destinations. The net effect is to reduce the number of 
functional cells. All of the above dependencies are hypothetical and illustrative 
only. In particular 






r(0(D(0/Co) 



(3.5) 



where in the latter ‘f{t) represents toxin input rate to the organ, and the 
remaining term v*(T(t)D(f)/Co)y|l + x'*T(f)j represents the toxin concentration 

rate-limited (Michaelis-Menton) reduction by functional cells. 

It has been observed that the realistic effect of toxin presence on A’, /i*, €t, etc. 

at time t may involve the entire past history of such exposure, e.g. 

A*(f)= A*(T(x),fl <x < f); the explicit function may be illustrated by the form 

A*(0 = exp^J^A*(T(x)<ix)j. A further model enhancement would recognize that 

toxin presence in the organ may result in cells that complete completing mitosis 



giving rise to damaged cells, e.g. possessing DNA adducts. These cells now 
become susceptible to repair. 
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3.2 Differential Equations For Deterministic Approximations 

Let Co denote the number of spaces/holes that cells may occupy. Then 

Co = Dit) + Git) + 2M(0 (3.6) 

since the cells in division potentially require two holes. Assume that if Cq » 1 
then 



D(f)/co ^ 0(0, G(o/co ^ r( 0 , m/co ^ m (3.7) 

in probability as Co ^ where the latter functions are 0(1). The following 
differential equations result from direct manipulation of (3.3): 

= -Ao(f ) - ea{t)y{t) + 2^(1 - o(f) - y(f ))^; (3.8) 

the last term occurs by virtue of (3.6). Also, 

^^ = Ao(f)-6o(f)y(f) (3.9) 

at 

^ = 5 = (3.10) 

dt 1 + kP 

We re-emphasize that all parameters are implicit functions of (current) toxin 
level, T(f), and possibly also absolute time, f, which can be viewed as the age of a 
mature organ, rather than an individual cell. This latter step is not taken in this 
paper. 

Note that in (3.8), (3.9), (3.10) the original parameters of (3.3) must be scaled: 
A*(T(0) is replaced by Xipit)), 6’(IX0) by 6ip(t)), /t*(T(f)) by //03(f)), v*(T(t)) by 
viPit)), T* by tCq, and k* by kCq. 

In Section 6 such scaling is exploited more extensively to deduce stochastic 
behavior of system state values. 
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3.3 Solution Without Toxin Input 

Suppose no toxin input exists, so T{t) = 0, Vf > 0. Then a steady-state solution 
to (3.8) and (3.9) may occur: set derivatives = 0 and solve the resulting equations 
to find 



The latter implies (i) that the larger the signal activity rate, 6, the smaller the 
ambient signal population, and the more quickly does signaling induce another 
cell to begin division; furthermore (ii) if 6 » X, the death rate, signifying highly 
efficient signaling, then a-^p.1 (2X + fi), which is equivalent to thinking of cells as 
behaving in pairs: once division completes the two daughters effectively compete 
to die (rate 2X). 

Although the equations (3.8) and (3.9) are non-linear and apparently cannot 
be solved explicitly one can obtain an approximate time-dependent solution as 
follows: 

(a) assume )(f) quickly reaches steady state, so 



This is a standard and useful approximation technique often invoked in 
biochemical kinetics problems that is variously abbreviated the quasi-steady-state 
assumption (QSSA) or the pseudo-steady-state hypothesis (PSSH); see Segel and 
Slemrod (1989) for a careful exposition of its rationale. The C^SA approach gives 
rise to the classical Michaelis-Menton used widely in pharmaceutical kinetics and 
compartment models; it already appears implicitly in our (3.5) and (3.10). 




(3.11) 



^'_ p{l-Xld) 



2X + P 



(3.12) 



y{t) = A/0; 0 < t 



(3.13) 
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(b) introduce (3.13) into (3.8) to find the approximation a{t): 



da{t)/dt = -(2A + ^)a + - X/d), 

the solution of which is 

a{t) = + [n(l - Xie)l{lX + /i)](l - j 

— > A/0)/(2A + /x) as f — 

as in (3.12). 



(3.14) 



(3.15) 



3.4 Dose-Response for Small Steady Toxin Exposure 

Of interest in risk analysis is the behavior of the dose response curve for small 
values of (toxic) dose. We approximate this by finding expressions for 



da{r) 



dr 



dy{T) 



T = 0' dr 



A 

and -9^ 
r = 0 dr 



r = 0' 



Suppose rit) = r = r ~ 0, a constant, and that exjx)sure has proceeded for some 

time so that a steady-state condition has been reached; this is modeled by letting 

da{t) _ dyjt) _ ^ _ q ( 3.8)- (3.10). To study the dose-response for t=0, 
dt dt dt 

differentiate (3.10) with respect to r using the function S in (3.5) (drop subscript 
"o" for convenience) 



0 



_ ^ dv dp aP vda P 



va dp ^ vap dp 



(3.16) 



dp dr I + kP dr 1 + kP 1 + kP dr (1 + ^p)^ dr 

Since )3 = 0 for t = 0, equation (3.6) yields for the (scaled) rate of toxin 
concentration at low exposure 



B 



T = 0 va 



1 XX •¥ n ^ t 

= 7 T\>0 

T = o e 



(3.17) 



This shows explicitly how low-dose toxin level in the organ tissue increases as a 
function of cell-cycle parameters X and /i, signaling efficiency, 6, and toxin 
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metabolic or binding rate v. All of these parameters eire evaluated at very low 
(zero) toxin levels. 

Next consider steady-state rate of change of scaled functional cell fraction, 
^ % = aa ■ this is the slope of a dose-(functioning cell count) response 

curve at very low levels of toxin. If (3.12) is differentiated and some 
rearrangements made, andif y= A/0as in (3.11), then 

“0 = [-^^ ^(2 + V/' + 2(1 - y )) + 00 (2A(1 - y)) + e^(2A + iiXnie)r] ^ ; . O.is) 

(2A + ^i) 

Since Px > 0, the sign of Ot is determined by the above linear combination of the 
derivatives of the cell cycle and signaling parameters A, fi, and 6 at t , and hence 
p, equal to zero. Note that if the illustrative relations (3.4) are invoked as given, 
then ap < 0, hence < 0 so a small increase in toxin induces a decline in the 
number of functioning cells; this is intuitively acceptable. However, a change in 
sign of any of the parameters' derivatives is capable of reversing the sign of the 
slope of the dose response curve, thus representing a tendency for hormesis at 
low dose. 

4. Models with Cell Aging 

In this section we generalize the models of the last section to allow cells to 
age, and hence exist in possibly different stages of fimctional effectiveness. 

We restrict attention to just two age stages, which we term new in) and old (o), 
distinguishing them by subscript. An arbitrary number of such age states can be 
introduced, but at the cost of increasing the number of parameters. 

Let Dn(t) (respectively Doit)) represent the number of new (respectively old) 
cells at time t. Mit) represents the number of dividing (mitotic) cells at time t and 
Git) the number of ambient signals or ghosts as before. The spatial constraint is 
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Co = Dnit) + Doit) + Git) + 2Mit). 



(4.1) 



Here are the transition probabilities assumed: 



Signaling Model with Age Dependence 
(Markov Generator) 



t t+dt 

Dnit), Do{t), ^ Dnit) - 1, Doit), Git) + 1,M(0 
Git), Mit); Tit) + 5*dt + odWit) 

Tit) (new cell death, signal creation) 

Dnit),Doit)-l,Git) + l,Mit) 

Tit) + S*dt + odWit) 

(old cell death, signal creation)! 

Dnit) - 1, Doit), Git) - l,Mit) + 1 
Tit) + 5*dt + odWit) 

(new cell signaled to divide, 
signal disappears) 

Dnit), Doit) - 1, Git) - 1, Mit) + 1 
Tit) + 5*dt + odWit) 

(old cell signaled to divide, 
signal disappears) 

Dnit) + 2, Doit), Git), Mit) - 1 
Tit) + 5*dt + odWit) 

(two newly-divided new cells 
emerge) 

^ Dnit) -1, Doit) + 1, Git), Mit) 
Tit) + 5*dt + odWit) 

(a new cell becomes an old cell) 



Probability 

XjiDniOdt 



AoDo(f)df 



6nGit)iDnit)/Co)dt 



elGit)iDoit)/Co)dt (4.2) 



H*Mit) dt 



<t>*Dnit)dt 



Once again all rates of cell transition, namely Xq, ol, p are presumed 
to be influenced by current toxin concentration. In turn, the latter is influenced by 
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the metabolic capability of the cells to remove toxin, 5 \T(t), Unit), Doit), tit)). 
When specific functions are required, we shall use for illustration 



s’im, D„(/),D„(o, <()) = T (<) - v; 


(4.3a) 


/ 

\ 




, P + ^^ = 1, P>0 


(4.3b) 






(4.3c) 




/(r(0)= A) 


(4.3d) 






(4.3e) 



The form of (4.3b) permits the adjustment of organ toxin concentration-caused 
cell death rate to be flexibly adjusted; if | < 0 the inner bracket must have its 

sign reversed: ( )■•■. As before, the above forms are hypothetical. 

4.1 Differential Equations for Deterministic Approximations 

Let Co denote the number of spaces/holes that cells may occupy. Then 



Co = Unit) + Doit) + Git) + 2Mit). 



Assume that 



Co Co Co Co 



(4.4) 



in probability as Co The following differential equations result from direct 
manipulation of (4.2) with scaled rate functions 

= -[K + ^n«n(0y(0 + "«n(0 " "oCO" /(O] (4.5) 

= -KcCoii) + ^«n(0- ^o«o(0y(0 (4.6) 
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= ~{d„an{t) + 0(,a(,{t)]y{t) + ^an{t) + 



(4.7) 



= a (t) — V a (t) — 



(4.8) 



These equations can be solved numerically, but not in closed parametric form. 



4.2 Solution Without Toxin Input 

Suppose no toxin input exists so T(t) = 0, t ^ 0. Then a steady-state solution to 
(4.5) - (4.7) may occur; set derivatives equal to 0 obtaining the following 
equations. 

0 = -(A„ + b)«n - 0„a„r + /i[l -a„-ao~Y] (4.9) 



0 = -^o<^o + 4>an~^o(^oY 


(4.10) 


0 = [— 2-o<^o- 


(4.11) 



Equation (4.10) yields 



= (y<P)[K + (4.12) 

Substituting (4.12) into equation (4.11) yields 

o=[-^«(V<^Po+M-^o]y+(^[(Vrt^ +M]+^o) (4.13) 

after dividing through by Oq- 
Rewriting (4.13) we obtain 



0= (-e„e„/0)r^ -y[e„(i/«K +e„ -A„(i/#)e„]+[(A„/#)+i]A„. (4.u) 

Thus Y satisfies a quadratic equation. Let Yo be the positive solution to the 
equation. If An(0) = 0, as in (4.3b) 




(4.15) 
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If Yo < 1, then y=Yo- If then ambient signals and not functional cells 

dominate the organ, which in reality would be long defunct even without toxin. 



rr 1 j 

If “T < 1 and 



< 


[14^1 


/ 


1 

1 

1 




'Pi 


\J 


L 



^0 

then Yo < 1. 

To find Oo equation (4.9) is used. Putci = (l/^)[Ao + Then 

~ ro] = [{Ki + + 0„)Ci + /i[ci + l]}«o = [iXo + /i[ci + 1] + 2A„ci}ao- 



(4.16) 



Thus, 



and 






(4.17) 



CCji — Cl OCq- 



(4.18) 



4.3 Dose-Response For Small Steady Toxin Input 

As mentioned earlier, the behavior of an organ-level dose-response curve for 
small values of toxic dose is of interest in risk analysis. In this section we indicate 
that the derivative of various cellular-level responses with respect to t may be 
evaluated at r= 0. 

If T = 0 then P = 0. Set the left-hand side of (4.8) equal to zero, then 
differentiate the right side with respect to t. After setting )3 = 0, this yields 

0 = 1 — (4.19) 



Thus 



a. 

dr 



r = 0 v„a„+Vo«o 



> 0 , 



as is physically plausible. 



(4.20) 
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It may be shown by further differentiation that under the specific 
parameterization of (3.4), or suitable alternatives, all derivatives of long-run 
cellular-level response can be explicitly evaluated at t = 0, as was shown earlier. 
The complex details are omitted here. 

5. A Model with Cell Aging and Organ Tissue Growth 

In this section we generalize the previous model to allow cells to 
spontaneously enter mitosis and multiply, thus effectively increasing the organ 
size, where the latter is defined as the number of cells in existence, in either state. 
In the previous section cell replication only occurred to replace other cells that 
previously have died, possibly as a result of toxin action. 

Consider the region to be occupied by cells to consist of spaces /holes that 
may be filled by cells. Let Co be the maximum number of such holes, in this case 
the number in an essentially mature organ. Initially not all such holes are active, 
for the organ is immature and hence growing. Hence let C(0 be the number of 
active holes at time t, while Cq - C(t) are currently quiescent or inactive. Now, 
generalizing the previous setup, let the current number of active holes /spaces be 

at) = Unit) + Doit) + Git) + 2M(0 (5.1) 

and augment the previous model to allow at) to gradually grow, although 
remaining bounded by Cq. 

5.1 Organ Growth Stimulated by Signaling 

In this model the presence of inactive holes /spaces effectively encourages the 
active cells, D„it) and Doit), to increase the number of signals and hence the 
number of active spaces capable of accomodating cells. We are implicitly 
modeling the competing effects of positive and negative, or inhibitory, growth 
factors. 
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The model is specified by the following transitions. 



Signaling-Driven Organ Growth 
(Markov Generator) 



t 



t+dt 



Probability 



Dn(t), Doit), D„(t), Do(t), Git) + 1, Mit), C(t) + 1 
Git), Mit), Tit) +6dt + odWit) 

at) Tit) (a quiescent cell space becomes 

active, signal creation) 

D„it) - 1, Doit), Git) + 1, M(f), at) 
Tit) +6dt + odWit) 

(new cell death, signal creation) 

^ Dnit), Doit) - 1, Git) + 1, M(f), at) 
Tit) +Sdt + odWit) 

(old cell death, signal creation) 

^ Dnit) + 2, Doit), Git) + 1, Mit) - 1, C(f) 
Tit) +Sdt + OdWit) 

(two newly-divided cells emerge) 

^ Dnit) - 1, Doit) + 1, G(f), Mit), at) 
Tit) + 6dt + odWit) 

(new cell becomes old) 






1 - 



C{t) 



dt 



XjjDnit)dt 



XoDoiOdt 



pi*Mit)dt 



(p^DniOdt 



'0 . 



(5.2) 



Dnit) - 1, Doit), Git) - 1, Mit) + 1, at) 
Tit) +6dt + odWit) 

(new cell signaled to divide, 
signal vanishes) 

Dnit), Doit) - 1, Git) - 1, Mit) + 1, at) 
Tit) +Sdt + odWit) 

(old cell signaled to divide, 
signal vanishes) 



0*Git)^^dt 
" ^ ’ C{t) 



OlGit) 



Ml 

C{t) 



dt 



18 



Once again, all rates of cell transition, namely, A^, 6^, n*, are 

presumed to be influenced by current toxin concentration. When specific 
functions are required we shall use 

(5.3) 

in addition to those of (4.3 a-e). 



Once again assume that as Cq 

Lq Cq Lq Cq 

in probability. 



(5.4) 

^0 



The following differential equations result from direct manipulation of (5.2). 

= -{h + 0]«n(O - + K'nit) - CCn{i) ~ «<j(0 ” HO) (5.5) 



dt I " "J "' ' n{t) 

= --^o«o (0 - ^0 y (0 + (0 



dt 

M0=_ 

dt 



ri{t) 

e„a„{t) + eoao{t) 



(5.6) 



??(0 



= [^n«n(0 + ^o«o(0][l - HO] 

^HO _ ^ 



dt 



y{t) + A„a„(f) + Ao«o(0 + [^o«<j(0 + ^«««(0P ” HO] (5-7) 

(5.8) 

(5.9) 



6. Stochastic Differential Equation Models 

The size of the state space makes the time-dependent behavior of the 
continuous-time Markov chain models of Sections 3-5 difficult to study. One 
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approach to studying the behavior is by Monte Carlo simulation; cf. Bois, et al. 
(1992). Another approach, adopted here, is to approximate the continuous-time 
Markov chain model by a diffusion process. Note that the continuous-time 
Markov chain model has absorbing states; for example, any state with D„(f) + 
Doit) = 0 is absorbing. Of course no self-respecting organ would ever start life in 
such a state, and a living organ would presumably die from other unmodeled 
causes long before such a state is reached, i.e. when D„it) + Doit) ^ D, some lower 
limit. Barbour (1976) discusses how long and over what ranges the vmderlying 
continuous-time Markov chain is approximated by a diffusion process of the type 
we derive; see also discussion to McNeil and Schach (1973). Stochastic 
differential equation models can be written for all the models of Sections 3-5. 
We will illustrate our approach by writing down a system of stochastic 
differential equations for the multivariate process (D„(f), Doit), Git), Cit), Tit)) for 
Markov generator (5.2). The system of stochastic differential equations is as 
follows; as explained subsequently, the dW-terms are normally distributed. 



dD„(() = -a;(t(())d„(i)a - e„{ni)y;{t){D„(i)/c(i))dt 

-i(T(l)}D„(t)dl + fl'(T(l)XC(l)-D„(t)-Dc(t)-G(l))dl 



^drifl =conditional^ 
expected change 
in (t,t + dt) 



+^e4p'(r(l)XC0)-D„(l}-DM-G(l))l,JW^(t] 



f diffusion = random^ 
"(^change in (t,t + dt) J 



(6.1) 
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(drift) 



dD„(l) = -A;(T{())D„(()rf(- e;(r(())G(()(D„(t)/C{t)>i<' 
+i(r(i))D„(i)di 

-^ee;{nt))G{t){D,(tyc{i)) (0 



►(diffusion) 



( 6 . 2 ) 



rfG(() = -e;(T(())G(f)(D„(()/C(())*-e;(r(t))G(tXD„(()/C(t)>i( 
+A;(r(/))D„(f)rf( + X'„(T{t))D„(t)dt 



+i;(no)D„(f) 



1-£W 

. Co J 



+^(r(())D„(o 



i-CW 

. Co J 



(drift) 



-^ee'„{T{l)}G(l){D„(l)/C(t))dWeJI) 



-^e0;(T(l))G(lXD,(l)/C(l)}dW^Jl) 
+^eX‘,(T(l))D,(t) rf (1) + J (() 






1 - 



C(f) 



i 



eC(T{‘)Po{l) 



1 - 



oo 

Co J 



dw^^it) 



dw^^{t) 



(diffusion) 



(6.3) 
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^ = «;(r(i))D„(f) 



,_£W 

. Co. 



i;(r(f))D„(() 



i-£W 



'0 J 



(drift) 



+j£^;(T(f))D,(0 



i_£(0 

. Q. 



dW^it) 



+j£^;(r(0)D„(0 



1_£(0 






(6.4) 



(diffusion) 



dTjt) 

dt 



<5'(r((),D„((),D„((),T'(f))ii 




+CT;(r((),D„((),D„((),i(0)rfw<,^(() 

(0 



{diffusion) 



(6.5) 



where {Wj„(()}-{Wi„(0}.{‘Ve,(0},{We^(()}.{%W}.{’^«.(')}'{'^{,WM’^<’r(')}. 
and {Wt(0) are independent standard Brownian motions; each dWit)-term is 
thus Gaussian with mean zero and variance dt. 

Setting the constant e = 1 lets the variances of the change in the number of 
active cells, ghosts in identified states, and active spaces be equal to the mean 
respective changes; this represents Poisson variability. Setting the constant £ > 1 
implies that the variability is larger than Poisson variability. This additional 
variability may be the result of inhomogeneities in the organ that are not 
explicitly modeled; cf. Bass, Robinson and Bracken (1978) for discussion of a 
distributed liver model with random inhomogeneities. Such overdispersion is 
frequently encountered in practice, and may be accounted for parsimoniously as 
we do; see McCullagh and Nelder (1983). 

We assume that as Co 



Ml 

Co 






^0 Cq Cq 



v{t), 



no 

Co 



Pit) (6.6) 
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in probability and rewrite D„(t), Doit), Git), Cit), and Tit) as follows. 

D„{() = Coa„(() + V^X„(l) (6.7a) 

D„(0 = Coa<,(() + VC^X„(() (6.7b) 

' G{l} = Cor{t] + ^Xg{t) (6.7c) 

C(t) = Co7)(f) + VQX^(() (6.7d) 

T{t) = Q,P{t) + ^Y{t) (6.7e) 

di (/) = Cot(/) + ^a^dW^it) (6.7f) 

The transition rates will be scaled and expanded as follows; 

^n(T(l)) = Kim/Co) = Kimh KWW>)/^I^ + 0(1) (6.8a) 

>‘l{n‘))=K{no/co)=Kim)+K{m)ni)/^+o{i) (6.8w 

m' (r(f)) = M(r(f)/Co ) = M V(<)) V(1)/ VQ + 0(1) (6.8c) 

enim) = e„(T(()/cb) = e„(i3{i))+K{m)ni)/^+om (6.8d) 

K(ni)) = eo(nO/Co) = eo(P(l)) + + 0(1) (6.8e) 

C(m)= um/co)=um)*i'n(m)n>)/^ +o(i) (6.8o 

lo'(r(<)) = &(r(0/Co) = um) + ^{m)y{‘)/yl^ + 0(1) (6.8g) 
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5’(D„((),D„((),r(0,r'(()) = Co5(D„(()/Co,D„(0/Co,r(()/Co-T'(t)/Co) 

= Co5(a„{t),aoit),p{t),t{t)) 



(6.8h) 



ir ^n(0 I ^ X' ^o(0 I ^ c' ^(0 I 

<’f(D»(<).D„((),r(().t'(f)) = V^<^7(D„(f)/Co,Do(')/Co-T(()/Cb,r(f)/Cb) 

(6.8i) 

= (a„ {t), Ooit), Pit), t ( 0 ) + 0 ( 1 ) 

Dividing equations (6.1) - (6.5) by Cq, scaling and substituting (6.7a) - (6.7f) 
results in a system of equations involving the deterministic components a„(0/ 
(Xoit), j(t), Tjit), and (Xt) and the corresponding stochastic components X„(f), XqH), 
Xgit), Xcit), and Y(t). See Appendix 1 for details. 

6.1 Second Moments 

The stochastic differential equations (A.l)- (A.IO) of Appendix 1, derived for 
the stochastic disturbances, can be written in matrix-vector form as 

dZit) = Mt)Zit) + Bit) dmt) (6.9) 

where in the present example Z is a 5x1 column vector, A is a 5x5 matrix, ^fW(t) is 
10x1 column of independent Gaussian noise terms and B(f) is a 5x10 matrix. 

It follows from theorem (8.5.5) of Arnold (1974) that 

^e[z?(()]= 2XA.y(f)E[z.(f)Zy(l)] + X(B.y(')f (SIO) 

j i 

and 
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j-iZi{t)Zj{t)] = XAit(OE[Zt(')Z;(')] 

k 

+IV(()E[Zt(()Z,(()] (6.11) 

jt 

k 

Thus, the approximate distribution of Ai(0, (respectively Doit), Git), Cit),Tit)) is 
normal with mean Coanit), (respectively CoOoit), Cojit) CoTlit), Copit)) and variance 
CoE[xJ(f)], tespectively Coe[x|(()],CoE|x|(()],CoE|x|(()],CoE|y 2(()]). The 
approximate covariance of Unit) and Doit) is CoE[Xnit)Xoit)], etc. 

We point out that exactly equivalent results can be obtained by introducing 
transforms for the state variables, computed by taking conditional mathematical 
expectations of exponential functions of the Markov generators. This is the route 
followed by McNeil and Schach (1973); see also Carpenter, Gaver, and Jacobs 
(1993). The transform approach provides a verification that the asymptotics 
described here actually lead to limiting Gaussian distributions and to Ornstein- 
Uhlenbeck processes. See Ethier and Kurtz (1986) for a more rigorous discussion. 
The approach taken here, while heuristic, is intuitively appealing and reaches the 
same ultimate conclusion. 

6.2 Numerical Examples 

In this section we present three numerical examples, 
fl. Constant Rate of Input of Toxin 

We assume chronic dosage of the organ by toxin that is delivered there at 
constant mean rate but with substantial randomness around that mean. 

Figures 1-5 depict the time variation of the expected numbers of new cells, 
of old cells, of total number of active (new plus old) cells, number of signals 
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(ghosts), and amount of toxin plus/minus two population standard deviations. 
There is no organ growth. The initial conditions are the limiting moments with 
no toxin input. The toxin input rate 0.2 units/ cell for all time. The new cells 
without toxin have mean time 100 before transitioning to old. The old cells have 
mean lifetime of 100 without toxin. The number of cells in the organ is 1.5x10^. 
The organ size is about that of a mouse liver. Since the cells in a mouse liver have 
a mean lifetime of about 200 days; cf. Bois, et al. (1992), the above is consistent 
with mouse liver behavior. The MATLAB 4*h and S**' order Runge-Kutta- 
Fehlberg numerical integration method was used to obtain the solutions. The 
MATLAB plotting algorithms were used to produce the figures. 

There is extra Poisson variability: e = 10^. 

Note that the total mean number of active (new and old) cells decreases as the 
result of toxin input, as seems intuitively reasonable. The mean number of new 
cells initially decreases, but eventually increases to a new steady-state value 
which is larger than the value with no toxin. The mean number of old cells 
initially increases briefly, but then decreases to a new steady-state value below 
the value for no toxin. The mean amount of toxin and the mean number of 
signals (ghosts) both increase to new steady-state values. 

The above hypothetical example illustrates how model parametric inputs 
translate into an account of the dynamics of subpopulations of cells. The details 
of the transition from one steady-state situation to another are of possible interest 
in that the organ may be in jeopardy during that transitional period because of, 
for example, new cell initial downward fluctuation. The effect on organ mortality 
is, however, not explicitly modeled. 
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b. A Toxin Pulse or Bolus Dose 

Figures 6-10 display the mean number of new cells, old cells, active cells 
(new and old), signals (ghosts) and amoimt of toxin plus and minus 2 standard 
deviations for the case in which the rate of toxin input is 0.2 per cell for the first 
50 time units and 0 thereafter. Once again there is no organ growth. The initial 
conditions are steady-state values when there is no toxin input. Note that the 
number of active cells decreases during toxin input, then increases and slightly 
overshoots the steady-state value with no toxin before returning to its steady- 
state value with no toxin. 

c. Organ Growth 

The model of Section 5.1 and (6.1) - (6.5) is initialized with the steady-state 
values for the second moments with no input of toxin. However, only 1/2 of the 
available spaces are active (may contain cells); the steady-state mean of new and 
old cells and ghosts for no toxin input are multiplied by 1 /2 and used as initial 
values. Figures 11-15 display the mean numbers of new, old, and active cells, 
ghosts, and active spaces plus/minus two standard deviations for the case of no 
toxin input and parameters ^n,0 = 0-5/ = 0.5, = 1/ ^o,l = 0.5. The mean 

number of active spaces is the maximum number by time 20. However, the mean 
number of old and new cells is still adjusting by time 100. The mean number of 
old cells initially decreases before increasing. The mean number of new cells 
initially increases before decreasing. 

The present model might be useful for describing the aftermath of a partial 
hepatechtomy: about one-half of a mature liver remains, and grows back with 
subpopulations of cells responding as shown. Note that, according to the present 
model in Figure 12, the organ must withstand an early signal to old cells to 
undergo substantial depletion to enter mitosis and bring forth replacements; the 
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latter are the new cells depicted in Figure 11. It may be that the sudden depletion 
of old cells to their low point at about f = 6 (Figure 12) actually puts the organ in 
jeopardy. 

Figures 16-21 display the mean number of new, old, and active cells, ghosts 
and active spaces plus /minus two standard deviations for the same parameters 
as used in Figures 11-15 but with two levels of toxin input. The dashed lines 
correspond to no toxin input. The solid lines correspond to a constant toxin input 
of 0.2 per cell per unit time. The graphs with positive toxin input have the same 
general shape as those with no toxin input. The mean number of new cells is 
initially less for a positive toxin input than for no toxin input. However, by time 
100 the mean number of new cells with positive toxin input is becoming larger 
than that for no toxin input. 

The effect of positive toxin input on the mean number of old cells is that it is 
initially larger than that for no toxin but by time 100 the mean number of old 
cells with positive toxin is less than that for no toxin. The effect of positive toxin 
input on the mean number of active (old and new) cells is to decrease it and 
delay its approach to a steady state value. The mean number of ghosts is larger 
with positive toxin input. Positive toxin input delays the overall growth in the 
number of active spaces. 
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mean number of new cells 



7 mean number of new cells plus minus 2 standard deviations 

x10 




Figure 1 
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mean number of old cells 



^ 7 mean number of old cells plus minus 2 standard deviations 
X 10 




Figure 2 
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mean number of active cells 



mean number of active ceils plus minus 2 standard deviations 
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Figiure 3 
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mean number of ghosts 



mean number of ghosts plus minus 2 standard deviations 



6 

X 10 




Figure 4 
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mean amt of toxin plus minus 2 standard deviations 




Figure 5 



mean number of new cells 



7 mean number of new cells plus minus 2 standard deviations 

X 10 




Figxire 6 



mean number of old cells 



7 mean number of old cells plus minus 2 standard deviations 

X 10 




Figure 7 
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mean number of active cells 



mean number of active cells plus minus 2 standard deviations 
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Figure 8 
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mean number of ghosts 



6 mean number of ghosts plus minus 2 standard deviations 

X 10 




Figxire 9 
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7 mean amt of toxin plus minus 2 standard deviations 

X 10 




Figxure 10 



mean number of new cells 



0.5 of the spaces are active 

mean number of new cells plus minus 2 standard deviations 




Figvire 11 
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mean number of old cells 



7 0.5 of the spaces are active 

X 10 mean number ot old cells plus minus 2 standard deviations 




Figure 12 
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mean number of active cells 



0.5 of the spaces are active initially 
I® mean number of active cells plus minus 2 standard deviations 




Figxire 13 
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mean number of active spaces 



0,5 of the spaces are active initially 

X 10 mean number of active spaces plus minus 2 standard deviations 




Figure 14 
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mean number of ghosts 



0.5 of the spaces are active initially 

mean number of ghosts plus minus 2 standard deviations 




Figure 15 
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mean number of new cells 



0.5 of the spaces are active initially 

mean number of new cells plus minus 2 standard deviations 




Figure 16 
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mean number of old cells 



0.5 of the spaces are active initially 

mean number of old cells plus minus 2 standard deviations 




Figure 17 
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mean number of old and new cells 
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Figure 18 
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mean number of ghosts 



0.5 of the spaces are active initially 
, 7 mean number of ghosts plus minus 2 standard deviations 

X 10 




Figiirc 19 
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mean number of active spaces 



0.5 of the spaces are active initially 

mean number of active spaces plus minus 2 standard deviations 




Figxire 20 



48 



16 



0.5 of the spaces are active initially 

mean amount of toxin plus minus 2 standard deviations 
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Figure 21 



A9 



Appendix 1 

This appendix provides the detailed asymptotic development of both the 
deterministic approximation (5.5) - (5.9) but also for the stochastic components 
introduced in (6.7). Start by introducing the normalization (6.7) into (6.1) - (6.5) 
and divide by Cq. The following results are obtained: 






dt 



^nW)) + 






(0 + (0 



dt 



0(1)) 

V^o_ 






dt 



[^(0-«n(0-«o(0-y(0] 



+;j^[x,(<) - x„(<) - x„{() - Xj(<)]U 






- On(0- «<.(<)- r('))| ‘IW^O) + 0 



VV'-u y 
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^o(0 + n^^^o(0 = ' 






dt 



Yff'iT 1 ^*^o(0+“7^^o(0 

^^oim^em))-^[rit^^ 

^(^( 0 ) + 









(A.2) 






. ^ ^s0M‘))r(i) “;dJ 4 vfey 
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dy{t) + 



1 




dXg{t) = 






V^o. 






rW+-i?-Xj{f) 






a„(f) + ^X„(f)' 

®o(0 ^o(0 



rff 









^oW'*' f^XQ{t) 

V^o 



«o(0+-i^Xo(f) rff 



um)*i'o{m)^ 

V'-o^ 

Wit))*KW‘))^ 

^ICo 

um)+K(m)^ 

V^o. 



--^^|£^n{P{t)]r{t)an{i)h{t)dWQ^{t^ 

V^o 



l-n(0-;^X,(f) 



rff 






«„(0 + - 7 ^ X „( f ) 









-;^T/^«o(/5(»))y(0“o(0/'i(0‘^'^«,(') 

+ (/5(f ))“»(') (t) + («())«„ (I) (() 

V^o 



+ V^&(«<)K(')[l - 1 (f)] “i W 4 , (() 



/ 

+ 0 

V 




(A3) 
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dn{t)-^^dXc{t) = 









dt 



^«(^( 0 )‘*‘^ n (/^( 0)^^^^(0 I ^ n ( 0 ’*"^^^^ n (0 I ^ ^(0 



dt 



(A.4) 



•gi^ V^^o(^(0)^o(0[i “ ^(0] (0 

/ 

V^^n(^(0)®«(0P “ ^(0] (0 + ^ 






+S'p(a„(t),aM-P{‘)M‘))^ 

+S-a^{a„(t).cco{t),P{t)M>))^ 

+^i„(“n(0'“o(0'W0'^(0) 






+ 




(T^rflV^(f) + 




The terms of order 1 in (A.l) - (A.5) yield the deterministic equations (5.5) - 
(5.9). 
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The terms of order 1/VQ in (A.l)- (A.5) give the following stochastic 
differential equations for the noise terms. 






-Sn(P{<M<) 



n(tr 



1 

v{t) 



Mmpcit)- X„ (0 - X,(f ) - X^ it)] 
(0 - «o (0 - y ( 0]^(0 






(A.6) 



h^£2^(/?(f ))(r?(t) -a„{t)-ao{t)-Y{t))dW^{t) 
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dXp(t) 

dt 



= -A„(/3(0)X„(() - A'(/3(())a„(()V(f) 
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-Uo(P{t))y{t)^dWe^{t) 
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dt 
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+A„(;8(())X„(() + Ai(«/))a„(()y{l) 
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(A.9) 
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